{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2452aa97",
   "metadata": {},
   "source": [
    "# Rotating domains"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dfab31b",
   "metadata": {},
   "source": [
    "We model configurations with rotating sub-domains, like electric motors, or wind mills. We generate independent meshes for the components, and glue them together using Nitsche's method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f544d9",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b22d004",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau = 0.0003\n",
    "tend = 1\n",
    "omega = 2*pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3de5928b",
   "metadata": {},
   "outputs": [],
   "source": [
    "square = MoveTo(0,0).Rectangle(1,1).Face()\n",
    "circo = Circle((0.5,0.5), 0.30001).Face()\n",
    "circ = Circle((0.5,0.5), 0.3).Face()\n",
    "holes = Circle((0.65,0.5), 0.05).Face() + Circle((0.35,0.5), 0.05).Face()\n",
    "\n",
    "square.edges.name=\"outer\"\n",
    "circ.edges.name=\"gammai\"\n",
    "holes.edges.name=\"hole\"\n",
    "circo.edges.name=\"gammao\"\n",
    "square.edges.name=\"wall\"\n",
    "square.edges.Min(X).name=\"inlet\"\n",
    "square.edges.Max(X).name=\"outlet\"\n",
    "outer = square-circo\n",
    "outer.faces.name = \"outer\"\n",
    "\n",
    "circ.faces.name = \"inner\"\n",
    "\n",
    "both = Compound([outer, circ-holes])\n",
    "mesh = Mesh(OCCGeometry(both, dim=2).GenerateMesh(maxh=0.05)).Curve(3)\n",
    "print (mesh.GetMaterials(), mesh.GetBoundaries())\n",
    "Draw (mesh);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8233e1fc",
   "metadata": {},
   "source": [
    "Define a GridFunction deformation describing the current configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7e50114",
   "metadata": {},
   "outputs": [],
   "source": [
    "fesdef = VectorH1(mesh, order=3)\n",
    "\n",
    "deformation = GridFunction(fesdef)\n",
    "defold = GridFunction(fesdef)\n",
    "\n",
    "def MeshRotation(angle, deform=None):\n",
    "    mesh.UnsetDeformation()\n",
    "    if not deform: deform = GridFunction(fesdef)\n",
    "    rotmat = CF( (cos(angle), -sin(angle), sin(angle), cos(angle))).Reshape( (2,2))\n",
    "    center = CF( (0.5, 0.5) )\n",
    "    pos = CF( (x,y) )\n",
    "\n",
    "    deform.Set( (rotmat-Id(2))*(pos-center), definedon=mesh.Materials(\"inner\"))\n",
    "    return deform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67ae9bde",
   "metadata": {},
   "outputs": [],
   "source": [
    "scene = Draw(mesh)\n",
    "tau1 = 1e-3\n",
    "for step in range(int(tend/tau1)):\n",
    "    MeshRotation(step*omega*tau1, deformation)\n",
    "    mesh.SetDeformation(deformation)\n",
    "    scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fad171c",
   "metadata": {},
   "source": [
    "Solve for a flow potential such that $\\frac{\\partial \\phi}{\\partial n}$ matches the normal velocity: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a6a3fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "fest = H1(mesh, order=3, dirichlet=\"inlet|outlet\")\n",
    "\n",
    "festgrad = VectorH1(mesh, order=3)\n",
    "gfutgrad = GridFunction(festgrad)\n",
    "\n",
    "ut,vt = fest.TnT()\n",
    "\n",
    "n = specialcf.normal(2)\n",
    "h = specialcf.mesh_size\n",
    "\n",
    "gfut = GridFunction(fest)\n",
    "    \n",
    "meshVelocity = (deformation-defold) / tau\n",
    "\n",
    "at = BilinearForm(grad(ut)*grad(vt)*dx)\n",
    "ft = LinearForm(fest)\n",
    "ft += -InnerProduct(meshVelocity,n)*vt*ds(definedon=\"hole\")\n",
    "\n",
    "contactt = ContactBoundary(mesh.Boundaries(\"gammai\"), mesh.Boundaries(\"gammao\"), volume=True)\n",
    "contactt.AddIntegrator (3/h*(ut-ut.Other())*(vt-vt.Other()))\n",
    "contactt.AddIntegrator (n*grad(ut)*(vt.Other()-vt)+n*grad(vt)*(ut.Other()-ut))\n",
    "\n",
    "\n",
    "def solveWind(gfut,at,ft):\n",
    "    contactt.Update (deformation, bf=at, intorder=10)\n",
    "\n",
    "    at.Assemble()\n",
    "    ft.Assemble()\n",
    "\n",
    "    # the solution field \n",
    "    gfut.Set((x), BND)\n",
    "    rt = ft.vec.CreateVector()\n",
    "    rt.data = ft.vec - at.mat * gfut.vec\n",
    "    gfut.vec.data += at.mat.Inverse(freedofs=fest.FreeDofs(), inverse = \"sparsecholesky\") * rt\n",
    "    gfutgrad.Set(Grad(gfut))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e465536c",
   "metadata": {},
   "outputs": [],
   "source": [
    "scene = Draw(gfut)\n",
    "tau1 = 2e-3\n",
    "for step in range(int(tend/tau1)):\n",
    "    defold.vec.data = deformation.vec\n",
    "    MeshRotation(step*omega*tau1, deformation)\n",
    "\n",
    "    mesh.SetDeformation(deformation)\n",
    "    solveWind(gfut,at,ft)\n",
    "    scene.Redraw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bee9d777",
   "metadata": {},
   "source": [
    "Solve a transport problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e9055a0",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fes = L2(mesh, order=3) \n",
    "u,v = fes.TnT()\n",
    "\n",
    "feshat = FacetFESpace(mesh, order=3)\n",
    "uhat, vhat = feshat.TnT()\n",
    "traceop = fes.TraceOperator(feshat, average=True)\n",
    "\n",
    "mesh.SetDeformation(MeshRotation(0))\n",
    "\n",
    "wind = -(meshVelocity - gfutgrad)\n",
    "\n",
    "a = BilinearForm(fes) \n",
    "a += -wind*u*grad(v)*dx\n",
    "uup = IfPos(wind*n, u, u.Other(bnd=0))\n",
    "a += wind*n*uup*v * dx(element_boundary=True) # upwind\n",
    "\n",
    "\n",
    "ahat = BilinearForm(feshat)\n",
    "\n",
    "f = LinearForm(fes)\n",
    "f.Assemble()\n",
    "\n",
    "gfu = GridFunction(fes)\n",
    "gfu.Set(exp(-10**2*((x-0.15)**2 +(y-0.5)**2)))\n",
    "\n",
    "solveWind(gfut,at,ft)\n",
    "scene = Draw(gfu, min=0, max=2, order=3, autoscale=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34c0f12a",
   "metadata": {},
   "outputs": [],
   "source": [
    "contact = ContactBoundary(mesh.Boundaries(\"gammao\"), mesh.Boundaries(\"gammai\"), volume=False)\n",
    "nc = contact.normal\n",
    "\n",
    "term = gfutgrad*nc * IfPos (gfutgrad*nc, uhat.Trace()*(-vhat.Trace().Other()), \\\n",
    "                               uhat.Trace().Other()*(vhat.Trace()))\n",
    "\n",
    "contact.AddIntegrator (term)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71247a73",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "t = 0\n",
    "cnt = 0\n",
    "deformation.vec[:] = 0\n",
    "w = gfu.vec.CreateVector()\n",
    "gfuhat = GridFunction(feshat)\n",
    "# what = gfuhat.vec.CreateVector()\n",
    "\n",
    "invm = fes.Mass(rho=1).Inverse()\n",
    "\n",
    "with TaskManager():\n",
    "    while t < tend:\n",
    "        defold.vec.data = deformation.vec\n",
    "        MeshRotation(t*omega, deformation)\n",
    "        \n",
    "        contact.Update (deformation, bf=ahat, intorder=10)\n",
    "        # apply the transport operator\n",
    "        mesh.SetDeformation(deformation)\n",
    "        solveWind(gfut,at,ft)\n",
    "        \n",
    "        gfuhat.vec[:] = traceop * gfu.vec\n",
    "        w[:] = a.Apply (gfu.vec) + traceop.T * ahat.Apply(gfuhat.vec)\n",
    "        \n",
    "        gfu.vec.data -= tau * invm * w\n",
    "        \n",
    "        if cnt%10 == 0:\n",
    "            mesh.SetDeformation(deformation)\n",
    "            scene.Redraw()   \n",
    "\n",
    "        t += tau\n",
    "        cnt +=1 \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b751277",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
